#####################
## run simulations ##
#####################
require(doSNOW)

source("../functions/pica2.R")
source("../functions/bounds_frechet_function_for_pica2.R")
source("../functions/pica1_use.R")

# code we use only for the simulation study 
source("bounds_frechet_function_new_estimand.R")
source("HPC_prepare.R")
source("HPC_loop_FINAL.R")

use_path <- "results/HPC_0724_big_"
ver_use  <- 1
SIM_SIZE <- 500 
testing  <- FALSE

## Design 
sample_size <- c(1500, 4500, 9000)
vio_list <- c(0, 0.05, 0.3)
choice_divergences <- c(rep(0, 3), seq(0, 1, 1/3))
outcome_divergences <- c(seq(0, 1, 1/3), rep(0, 3))
divergences <- seq(1:7)

design <- expand.grid(sample_size, vio_list, divergences)
colnames(design) <- c("sample_size", "vio", "divergence")

sim_total <- SIM_SIZE*nrow(design)

# Number of MC draws
n_dir_MC <- 100
n_norm_MC <- 10

## effects.mat 
J <- 3
effects.mat <- matrix(FALSE, J, J)
dimnames(effects.mat) <- list(ctrl = NULL, treat = NULL)
effects.mat[2, 1] <- TRUE
effects.mat[3, 1] <- TRUE
effects.mat[3, 2] <- TRUE

## #####################
## run on HPC Della
## ######################
# Read job array info
job_idx <- as.numeric(Sys.getenv("SLURM_ARRAY_TASK_ID"))
n_jobs  <- as.numeric(Sys.getenv("SLURM_ARRAY_TASK_COUNT"))

if(testing == TRUE){
  job_idx <- 1
  n_jobs  <- 1
}

# Timing
start_time <- Sys.time()
run_ts <- paste0(strftime(start_time, "%y%m%d"), "_ver_", ver_use)

# total simulation size
total_sim_size <- SIM_SIZE * nrow(design)

# Map from job_idx to sim_num
n_sims  <- total_sim_size
n_tasks <-  floor(n_sims/n_jobs)
remainder <- n_sims %% n_jobs
jobs <- (n_tasks*(job_idx-1)+1):(n_tasks*(job_idx))
if (job_idx<=remainder) {
  jobs <- c(jobs, n_jobs*n_tasks+job_idx)
}

if(testing == TRUE){
  jobs <- 1:20
}

# Run
for (idx in jobs){
  out <- tryCatch(
    {HPC_loop(s = idx)},
    error=function(e){
      message(paste0("ERROR: ", e))
      return(NULL)},
    warning=function(w){
      message(paste0("WARNING: ", w))
      return(NULL)}
  )
}